NONLINEAR  WAVE  INTERACTION  WITH  SUBMERGED  AND 
SURFACE-PIERCING  POROUS  STRUCTURES 


Okey  Nwogu1  and  Zeki  Demirbilek2 


A  coupled  Boussinesq-boundary  integral  method  is  developed  to  simulate  nonlinear  water 
wave  interaction  with  structures  consisting  of  multiple  layers  with  different  physical  and 
hydraulic  characteristics.  The  flow  field  in  the  water  region  is  modeled  with  a  modified 
set  of  Boussinesq-type  equations,  with  additional  terms  to  account  for  the  flow  of  water 
into/out  of  the  porous  region.  The  equations  of  motion  for  the  porous  regions  include  an 
empirical  Forchheimer-type  term  for  laminar  and  turbulent  frictional  losses,  and  an 
inertial  term  for  acceleration  effects.  A  boundary  integral  formulation  based  on  Green’s 
third  identity  is  used  to  close  the  problem  for  the  porous  region.  The  coupled  equations 
for  the  evolution  of  the  free  surface  and  boundary  values  of  the  tangential  velocities  are 
integrated  in  time  using  an  iterative  Crank-Nicolson  scheme.  At  each  time  step,  the 
Boussinesq  problem  is  solved  for  the  water  region  to  determine  the  pressure  at  porous 
interface.  The  boundary  integral  problem  for  the  porous  region  is  then  solved  to 
determine  the  normal  velocities  along  porous  interface.  The  model  is  used  to  investigate 
wave  interaction  with  a  vertical  surface-piercing  porous  structure  and  wave  transmission 
over  submerged  breakwaters.  Comparisons  between  the  numerical  model  predictions  and 
laboratory  data  show  generally  good  agreement  for  both  the  wave  field  inside  the 
structure  and  the  reflection/transmission  coefficients. 


1.  INTRODUCTION 

Porous  structures  such  as  rubblemound  breakwaters  typically  consist  of 
multiple  dissipative  layers  with  different  physical  and  hydraulic  characteristics. 
A  better  understanding  of  the  dissipation,  reflection  and  transmission 
characteristics  of  porous  structures  is  required  to  accurately  predict  wave 
conditions  in  navigation  channels  and  harbors,  sediment  dynamics  behind 
submerged  breakwaters,  and  erosion  and  sedimentation  processes  near  coastal 
inlets. 

Early  theoretical  studies  of  wave-porous  bed  interaction  were  concerned 
with  the  attenuation  of  surface  waves  propagating  over  permeable  horizontal 
seabeds  (e.g.  Putnam,  1949;  Reid  and  Kajiura,  1957).  Linear  wave  theory  was 
used  to  describe  the  overlying  wave  field  while  Darcy’s  law  was  used  to 
describe  the  flow  inside  the  porous  layer.  Darcy’s  law  is  valid  for  flows  with 
low  Reynolds  number.  As  the  flow  velocity  or  particle  size  increases,  Darcy’s 
linear  relationship  between  the  discharge  velocity  and  pressure  gradient  breaks 
down  and  the  quadratic  drag  law  proposed  by  Forchheimer  (1901)  becomes 
more  applicable.  Sollitt  and  Cross  (1972)  presented  a  modified  set  of  equations 
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for  fluid  motion  within  a  porous  medium  that  included  both  a  Forchheimer-type 
drag  term  and  inertial  effects.  An  eigenfunction  expansion  method  was  used  to 
solve  the  linearized  form  of  the  equations  for  waves  normally  incident  on 
vertical  surface -piercing  porous  structures. 

In  many  applications,  it  is  important  to  investigate  wave  propagation  over 
non-horizontal  seabeds  or  wave  interaction  with  non-rectangular  breakwaters. 
The  porous  medium  could  also  consist  of  multiple  layers  with  different 
porosities  and  dissipative  properties.  Several  investigators  (e.g.  Sulisz  1985; 
Mallayachari  and  Sundar  1994)  have  used  linear  boundary  integral  methods  to 
simulate  wave  interaction  with  multi-layer  porous  structures.  Linear  models, 
however,  cannot  accurately  describe  the  nonlinear  wave  steepening  and 
breaking  processes  that  occur  on  the  slopes  of  submerged  and  surface -piercing 
trapezoidal  structures.  These  complicated  processes  can  be  better  addressed  by 
numerically  solving  either  the  full  Navier-Stokes  equations  (e.g.  Liu  et  al.  1999) 
or  asymptotic  theories  for  weakly  nonlinear/weakly  dispersive  waves. 

Flaten  and  Rygg  (1991)  developed  an  asymptotic  Boussinesq-type  theory 
for  shallow  water  wave  propagation  over  a  permeable  bed.  Darcy’s  law  was 
used  to  model  the  flow  inside  the  porous  layer.  Cruz  et  al.  (1997)  developed  a 
two-layer  weakly  nonlinear  Boussinesq  model  that  included  both  fluid 
acceleration  effects  and  quadratic  drag  losses  inside  the  porous  layer.  The 
dispersive  and  damping  properties  of  the  equations  were  also  improved  by 
introducing  additional  dispersive  terms  to  the  momentum  equations  using  the 
approach  of  Madsen  and  Sorensen  (1992).  Flsiao  et  al.  (2002)  extended  the  two- 
layer  approach  to  large  amplitude  waves  by  using  a  fully  nonlinear  set  of 
Boussinesq  equations  for  the  water  region.  While  the  two-layer  Boussinesq 
approach  is  suitable  for  modeling  wave  propagation  over  permeable  seabeds,  it 
is  not  ideal  for  studying  wave  interaction  with  multi-layer  porous  structures  of 
arbitrary  cross-section.  The  equations  assume  that  the  thicknesses  of  the  porous 
layers  are  much  smaller  than  the  incident  wavelength  and  the  interfacial/seabed 
slopes  are  mildly  varying.  For  structures  that  are  a  fraction  of  a  wavelength,  the 
spatial  discretization  schemes  used  for  the  open  water  region  are  also  too  coarse 
to  resolve  details  of  the  flow  field  inside  the  porous  structure. 

In  this  paper,  a  numerical  model  is  presented  for  simulating  nonlinear  wave 
interaction  with  multi-layer  porous  structures  of  arbitrary  shape.  The  coupled 
model  consists  of  a  set  of  fully  nonlinear  Boussinesq  equations  for  the  water 
region  and  boundary  integral  equations  for  each  of  the  porous  layers.  This 
approach  does  not  make  any  assumption  as  to  the  profile  of  the  flow  field  inside 
the  porous  medium.  The  numerical  model  is  initially  calibrated  with  data  from 
laboratory  experiments  on  wave  interaction  with  a  vertical  porous  breakwater 
(Requejo  et  al.  2002).  The  numerical  model  is  then  used  to  investigate  wave 
transmission  over  submerged  breakwaters  and  compared  with  the  laboratory 
data  of  Seabrook  and  Flail  (1998). 
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2.  THEORETICAL  BACKGROUND 

Governing  Equations  for  Water  Region 

Consider  the  interaction  of  weakly  dispersive  surface  waves  with  a  multi¬ 
layered  surface-piercing  breakwater  as  shown  in  Figure  1.  A  fully  nonlinear  set 
of  Boussinesq  equations  is  used  to  describe  the  fluid  motion  in  the  water  region. 
The  evolution  of  the  free  surface  r](x,t)  is  governed  by  the  depth-integrated  mass 
conservation  equation: 

r)  ,+V-g  =  q(np)  (1) 

where  V  =  (d/dx ,d/dy),  r/(x,i)  is  the  water  surface  elevation,  (^represents  the 

normal  flow  through  the  water-porous  layer  interface  and  Q(x,t )  is  the  volume 
flux  density  given  by: 

Q  =  \n_u(w)  dz  (2) 

is  the  horizontal  fluid  velocity  in  the  water  region  and  h(x)  is  the  water 
depth. 


Figure  1.  Definition  Sketch. 

The  momentum  equation  is  derived  by  evaluating  the  Euler  equations  at  the  free 
surface  to  obtain: 

tif^l  +  wt(w)|  V//  +  gV  ij  +  (u^v)  ■  v)^’1'1  +  vvJw)VwJw)  =  0  (3) 

where  vv<M)  is  the  vertical  velocity.  The  system  of  equations  (1)  to  (3)  requires 
knowledge  of  the  vertical  profile  of  the  flow  field  in  order  to  evaluate  the 
horizontal  and  vertical  velocities  at  the  free  surface  and  the  volume  flux  density. 
For  shallow  water  waves,  the  vertical  profile  of  the  flow  field  is  obtained  by 
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expanding  the  velocity  potential,  <D,  as  a  Taylor  series  about  an  arbitrary 
elevation,  za,  in  the  water  column  (Nwogu  1993): 


<D  (v\x,z,t)  =  c/>a+(z-za) 


^  +^[(zff+/?)2-(z  +  /7)2]vVa(4) 


where  <f>  =  <b{"']  (x,z,t)  ■  The  vertical  velocity  at  the  bottom  of  the  water  region 

depends  on  the  inflow/outflow  from  the  porous  region  and  is  obtained  by 
enforcing  the  continuity  of  normal  velocity  at  the  water-porous  layer  interface: 


By  substituting  Eq.  (5)  into  Eq.  (4)  and  differentiating,  we  obtain  the  vertical 
distributions  of  the  horizontal  and  vertical  velocities  as: 


u(w\x,z,t)  =  =  ua+(za-z)[v(ua-Vh)  +  (y.ua)Vh-Vq^]  (fi) 

+  |[(z«+/02-(^)2k(v-«o) 


w(w\x,z,t)  =  - -  =  q[p)  -  \ua  ■  Vh  +  (z  +  h)V  •«  1 

where  the  velocity  variable  ua  =  u(x,za,t)  is  the  horizontal  velocity  at  z  =  za.  The 
volume  flux  density  required  to  evaluate  the  mass  conservation  equation  is 
obtained  by  substituting  Eq.  (6)  into  Eq.  (2)  and  integrating  to  yield: 

Q  =  (h  +  q) !*„+  (Zg^)  V(V-ng)+  (8) 

[v(ua  -Vh)  +  (V -ua)Vh- Vq[p) ]  J 

The  momentum  equation  can  be  written  in  terms  of  the  velocity  variable  ua  by 
substituting  the  expressions  for  uiw)  and  wM  (Eqs.  6  and  7)  into  Eq.  3  to  obtain: 

u„t  +  gVn  +[<>^k>  +  +  (za-4Mua,  •  VA)+ (V- ua,)S7h -q^\ 

+  -^a+h)2-(h+nfW-uJ  -  k, .M+ih+rfF-u^-^ 

+1]-fA”\<\+1]-V{v{h+rm-ua)  =  0 
h+rf  1  1  /z+77 


(9) 
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where  v  and  fw  are  the  eddy  viscosity  and  bottom  friction  coefficient  respectively 
intended  to  account  for  wave  energy  dissipation  due  to  wave  breaking  and 
bottom  friction.  The  wave  breaking  model  is  based  on  the  eddy  viscosity  model 
of  Nwogu  (1996). 

Equations  (1),  (8)  and  (9)  represent  a  modified  set  of  fully  nonlinear 
Boussinesq  equations  for  the  water  region.  Compared  to  the  fully  nonlinear 
equations  of  Nwogu  (1996),  Eqs.  (1)  and  (9)  contain  additional  terms  involving 
q(np)  to  account  for  flow  into  and  out  of  the  porous  structure.  The  normal  flow 

through  the  water-porous  layer  interface  will  be  determined  by  solving  a 
boundary  value  problem  for  the  porous  region. 


Governing  Equations  for  Porous  Region 

The  asymptotic  Taylor  series  expansion  used  to  describe  the  flow  field  in 
the  water  region  could  also  be  applied  to  shallow  porous  layers  (e.g.  Hsiao  et 
al.,  2002).  However,  this  approach  is  cumbersome  to  apply  to  structures  that 
consist  of  multiple  porous  layers  of  arbitrary  cross-sections.  Since  the  flow 
inside  the  porous  layer  satisfies  the  Laplace  equation,  Green’s  third  identity  may 
be  applied  over  any  closed  boundary  S  to  relate  the  values  of  the  potential  inside 
the  bounded  region  with  the  boundary  values  of  the  velocity  potential  ^'’((f) 
and  its  normal  derivative  50^’  (£)/dn.  For  field  points  x  that  lie  along  the 
boundary,  the  velocity  potential  is  given  by: 


<S(rt(x)  = 


5<5(p) 

dn 


(f)-®(rt(f)^(*;f) 

on 


dS 


(10) 


where  G  =  -In  |x-f|  is  the  Greens  function  and  c  represents  the  source  point  on 
the  boundary  S  over  which  the  integration  is  performed.  The  integral  equation  is 
evaluated  numerically  by  dividing  the  porous  layer  boundaries  into  a  finite 
number  of  segments  with  the  values  of  0(/I)  and  d<l),pl/dn  assumed  to  be  constant 
over  each  segment  and  defined  at  the  midpoints  of  the  segments.  The  integral 
equation  for  each  porous  region  can  thus  be  rewritten  as  a  matrix  equation  in 
terms  of  discrete  values  of  cE></J)  and  dA'>lp)ldn\ 


where  N  is  the  number  of  segments  on  the  boundary  and  the  matrix  coefficients 
Ay,  By  correspond  to  the  integral  of  the  Green’s  function  or  its  normal  derivative 
over  the  individual  segments.  When  i  ^  j,  the  integrals  are  approximated  by 
assuming  the  Green’s  function  (or  normal  derivative)  to  be  constant  over  each 
segment  and  evaluated  at  the  midpoint  of  the  segment.  When  i  =  j.  the  integrals 
become  singular  and  are  evaluated  analytically. 

Each  porous  layer  may  have  up  to  four  different  types  of  boundary 
conditions  corresponding  to:  1)  matching  boundary  with  water  region,  2)  free 
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surface,  3)  matching  boundary  with  another  porous  layer,  and  4)  impermeable 
boundary. 

Along  matching  boundaries  with  the  water  region,  we  enforce  the 
continuity  of  pressure  and  normal  velocities.  The  solution  of  the  problem  for  the 
water  region  provides  pressure  boundary  conditions  for  the  porous  layer  while 
the  solution  for  the  porous  layer  provides  normal  velocity  boundary  conditions 
for  the  water  region.  Following  Sollitt  and  Cross  (1972),  we  assume  that  the 
flow  in  the  porous  region  is  governed  by  a  Forchheimer-type  resistance  term 
plus  an  inertial  term  to  account  for  fluid  acceleration  around  solid  particles.  The 
continuity  of  pressure  boundary  condition  is  obtained  by  evaluating  the 
modified  Euler  equations  for  the  porous  region  along  the  interface  z  =  z;.  This 
can  be  expressed  in  terms  of  the  tangential  velocity  along  the  interface 
q\p>  =  Va>('’)  as: 

cmq®  +  gVz,.  +  y  («,W  •  V)i»  +  £=- 

+  fM?+fM(?Up>-<>  +  ^>2  =  --Mp?  (12) 

p 


where  c,„  is  an  inertial  coefficient,  which  can  be  written  in  terms  of  the  added 
mass  coefficient,  ca,  using  cm  =  1  +  (1-1)  ca,  1  is  the  effective  porosity  of  the 
medium,  (u\p\w\p))  are  the  discharge  velocities  on  the  porous  layer  side  of  the 


interface,  and  f  and  f  are  empirical  laminar  and  turbulent  friction  factors.  The 
friction  factors  can  be  related  to  the  permeability  of  the  medium  (Sollitt  and 
Cross,  1972)  or  the  characteristic  stone  size,  d,  and  porosity  by  (e.g.  Van  Gent 
1995): 


f, 


f, 


1  + 


7.5 

K 


c  j 


d 2 

|(1-A)  1 

T3  d 


(13) 


ao,  and  h{]  are  empirical  constants,  and  Kc  =  id H(/xl)  is  the  Keulegan-Carpenter 
number  with  H  the  characteristic  wave  height.  The  empirical  coefficients 
depend  on  the  particle  shape,  surface  roughness  and  Reynolds  number  of  the 
flow. 

The  pressure  boundary  condition  (Eq.  12)  is  also  applied  along  phreatic 
surface  boundaries  in  the  porous  region  with  zt  =  t]  and  pi  =  0.  The  location  of 
the  phreatic  surface  is  obtained  by  solving  the  kinematic  free  surface  boundary 
condition:  ,  . 

1  50(p)  r— — — 

"'=I^T'/1+V"'V"  <i4> 

The  normal  velocities  5O(p)/0«  are  obtained  as  part  of  the  solution  of  the 
boundary  integral  equation. 
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Along  matching  boundaries  between  porous  layers,  we  enforce  the 
continuity  of  <DW  and  oO^Vo'm  while  the  zero  normal  flow  boundary  condition 
=  0)  is  enforced  along  impermeable  boundaries. 

Numerical  Solution 

The  governing  mass  and  momentum  equations  for  the  evolution  of  the 
free/phreatic  surfaces  and  velocity  variables  are  integrated  in  time  using  a 
modified  Crank-Nicolson  scheme.  At  each  time  step,  we  initially  solve  the 
Boussinesq  problem  for  water  region  to  determine  the  pressure  forcing  along  all 
porous  interfacial  boundaries.  The  porous-region  momentum  equation  is  then 
integrated  to  determine  the  velocity  potential  along  the  interface  while  the 
kinematic  boundary  condition  is  integrated  to  update  the  location  of  the  phreatic 
surface.  The  boundary  integral  problem  for  the  porous  region  is  then  solved  to 
determine  normal  velocities  along  the  interfacial  boundaries.  Finally,  the 
Boussinesq  problem  is  re-solved  with  updated  values  of  the  normal  velocity 
along  porous  interface.  This  procedure  is  repeated  at  each  time  step  until 
convergence. 

3.  NUMERICAL  RESULTS 

Wave  Interaction  with  a  Vertical  Surface-Piercing  Breakwater 

The  performance  of  the  numerical  model  was  initially  evaluated  with  data 
from  the  laboratory  experiments  of  Requejo  et  al.  (2002),  who  investigated  the 
interaction  of  periodic  waves  with  a  crib-style  rectangular  breakwater.  The 
experiments  were  conducted  in  a  12-m  long  and  0.56-m  wide  wave  flume  using 
two  structure  widths  (B  =  0.56m,  1.12m),  two  characteristic  stone  sizes  (d  = 
0.0105m,  0.0378m),  three  wave  periods  (T  =  1.2s,  1.7s,  2.7s),  and  five  wave 
heights  (II  =  0.02m  to  0.09m).  Five  wave  probes  were  used  to  measure  the  wave 
field  in  front  of  the  structure  and  two  or  three  wave  probes  were  used  to 
measure  inside  the  structure.  The  experimental  layout  is  shown  in  Fig.  2. 
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Figure  2.  Layout  of  laboratory  experiments  of  Requejo  et  al.  (2002). 


Table  1.  Summary  of  Best-Fit  Empirical  Constants  for  H=0.034m,  7=1 .7s. 
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Figure  3.  Comparison  of  measured  and  predicted  wave  heights  for 
H=0.034m,  7=1. 7s. 


We  investigated  the  sensitivity  of  the  numerical  model  predictions  of  the  wave 
heights  inside  the  porous  structure  to  the  values  of  the  empirical  constants  a0,  ba. 
Numerical  simulations  were  carried  out  for  one  of  the  test  conditions  ( H  = 
0.034m,  T  =  1.7s)  for  different  combinations  of  a0  and  bm  with  a0  varying  from 
0  to  300,  and  ba  from  0  to  3.  These  range  of  values  were  selected  to  cover  the 
optimal  values  a0=  100  and  ba  =1.1  recommended  by  Requejo  et  al.  (2002). 
The  value  ca  =  0.4  was  used  for  added  mass  coefficient  as  suggested  by  Van 
Gent  (1995).  The  values  of  a0  and  ba  that  yielded  the  least  square  error  between 
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the  measured  and  predicted  wave  heights  inside  the  structure  are  summarized  in 
Table  1.  The  empirical  constant  b0  associated  with  turbulent  losses  was  smaller 
for  Structures  #3  and  #4  that  have  the  larger  stone  size. 

The  numerically  simulated  wave  height  distributions  with  the  best-fit 
coefficients  are  compared  to  the  measured  wave  heights  in  Figure  3.  Although 
the  reflection  coefficients  were  similar  for  some  of  the  test  conditions,  the 
different  structures  have  fairly  different  wave  height  variations  inside  the 
structure.  Excellent  agreement  is  observed  between  the  measured  and  predicted 
wave  heights  inside  the  structure  for  Structures  #1,  #2,  and  #3. 


o 


o 


H/gT 


H/gT 


Structure  #3 

— • —  Measured 

— A —  Predicted 

Structure  #4 


H/gf  H/gT 

Figure  4.  Measured  and  predicted  reflection  coefficients  (T  =  1.7s). 


The  best-fit  empirical  constants  for  Structures  #3  and  #4  (a0  =  100,  ba  =  0.6) 
were  used  in  numerical  simulations  for  other  wave  conditions  and  structure 
configurations  to  investigate  the  potential  relative  error  in  using  a  universal  set 
of  constants  to  estimate  the  reflection  properties  of  porous  structures.  The 
measured  and  computed  reflection  coefficients  are  shown  in  Figure  4  for  T 
=1.7s,  and  Figure  5  for  T  =  1.2s.  Excellent  agreement  is  observed  between  the 
predicted  and  measured  reflection  coefficients  for  Structures  #1,  #2  and  #3  for  T 
=  1.7s.  The  numerical  model  also  reproduces  the  observed  trends  in  the 
variation  of  reflection  coefficient  with  wave  steepness.  For  the  T  =  1.2s  test 
conditions,  the  numerical  model  predicts  larger  values  of  the  reflection 
coefficient  (5-25%)  but  generally  reproduces  the  correct  trends  in  the  variation 
of  reflection  coefficient  with  wave  steepness. 
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Figure  5.  Measured  and  predicted  reflection  coefficients  (7  =  1.2s). 


Wave  Propagation  over  a  Submerged  Breakwater 

Seabrook  and  Hall  (1998)  carried  out  an  extensive  series  of  tests  to 
investigate  wave  transmission  over  submerged  rubblemound  breakwaters.  The 
tests  were  conducted  in  a  47-m  long  and  1-m  wide  wave  flume  at  Queens 
University,  Canada.  The  layout  for  the  experiments  is  shown  in  Figure  6.  The 
breakwater  cross-section  had  a  crest  height  of  hs  =  0.25m  and  1V:3H  side 
slopes.  The  breakwater  consisted  of  a  relatively  coarse  inner  core  layer  with  dsn 
=  0.017m  and  two  layers  of  primary  armor  with  d  =  0.059m.  Over  800  irregular 
wave  tests  were  performed  with  13  breakwater  geometries,  5  water  levels,  3 
spectral  peak  periods,  and  4  significant  wave  heights. 

Since  the  transmission  coefficient  for  submerged  breakwaters  depends  on 
parameters  such  as  the  relative  submergence  depth  and  relative  crest  width,  we 
carried  out  two  sets  of  numerical  simulations.  The  first  set  of  simulations  were 
conducted  for  a  constant  crest  width  B  =  0.6  m  and  varying  submergence  depths 
ds  =  0.05,  0.1,  0.15  and  0.2m.  For  the  second  set  of  simulations,  we  varied  the 
crest  width  B  =  0.6,  1.5,  2.5  and  3.5m  and  used  a  submergence  depth  of  0.05m. 
The  numerical  simulations  were  performed  for  an  incident  sea  state 
characterized  by  a  JONSWAP  spectrum  with  Hmo  =  0.1  m,  Tp  =  1.5  s  and  y  = 
3.3.  A  porosity  of  0.5  and  empirical  constants  a0  =  100,  ba  =  0.6,  and  ca  =  0.4 
were  used  for  the  model  simulations. 
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Figure  6.  Layout  for  the  Seabrook  and  Hall  (1998)  experiments. 

Figure  7  shows  a  comparison  of  the  variation  of  the  measured  and  predicted 
transmission  coefficients  with  relative  submergence  depth,  ds.  For  non-breaking 
waves  with  ds/Hmo  =  2,  the  transmitted  wave  height  is  about  80%  of  the  incident 
wave  height.  As  ds  decreases,  the  transmitted  height  decreases  to  about  40%  of 
the  incident  wave  height.  Wave  breaking  over  the  breakwater  also  plays  an 
important  role  in  the  wave  energy  dissipation  process.  Figure  9  shows  a 
comparison  of  the  wave  transmission  coefficient  as  a  function  of  the  breakwater 
crest  width  for  a  submergence  depth  of  0.05m.  The  results  show  a  decrease  in 
transmission  coefficient  with  increasing  crest  width  due  to  both  wave  breaking 
and  energy  dissipation  inside  the  porous  structure. 


Figure  7.  Variation  of  transmission  coefficient  with  relative  submergence 

depth. 
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Figure  8.  Variation  of  transmission  coefficient  with  structure  crest  width. 
CONCLUDING  REMARKS 

A  coupled  Boussinesq-boundary  integral  model  has  been  developed  to 
simulate  nonlinear  wave  interaction  with  surface-piercing  and  submerged 
porous  structures.  The  model  solves  vertically  integrated  Boussinesq-type 
equations  for  the  water  region  and  evolution  equations  for  the  phreatic  surface 
and  boundary  values  of  the  tangential  velocities  in  the  porous  region.  A 
boundary  integral  method  is  used  to  close  the  problem  for  porous  region  and 
determine  the  relationship  between  the  normal  and  tangential  velocities  along 
the  boundaries.  Since  the  energy  dissipation  within  porous  structures  is 
governed  by  empirical  laminar  and  turbulent  frictional  factors,  the  applicability 
of  the  model  formulation  and  determination  of  the  friction  factors  have  been 
investigated  by  comparing  numerical  model  predictions  of  wave  reflection  and 
transmission  coefficients  through  porous  breakwaters  with  data  from  two 
laboratory  experiments.  Overall,  the  model  performed  reasonably  well  and  was 
able  to  reproduce  the  observed  variation  of  reflection  coefficient  with  wave 
steepness  with  one  set  of  empirical  constants.  Good  agreement  was  also 
observed  for  the  wave  field  inside  the  structure,  provided  the  empirical 
coefficients  are  suitably  chosen. 
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